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We calculate the nonequilibrium local density of states on a vibrational quantum dot coupled to 
two electrodes at T = 0 using a numerically exact diagrammatic Monte Carlo method. Our focus 
is on the interplay between the electron-phonon interaction strength and the bias voltage. We find 
that the spectral density exhibits a significant voltage dependence if the voltage window includes one 
or more phonon sidebands. A comparison with well-established approximate approaches indicates 
that this effect could be attributed to the nonequilibrium distribution of the phonons. Moreover, 
we discuss the long transient dynamics caused by the electron-phonon coupling. 
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Recent experiments in the field of molecular electronics 
have pointed out the importance of electron-phonon in¬ 
teractions for the charge transport on the nanoscale»i - — 
In these experiments, a nanostructure - such as a sin¬ 
gle molecule or a carbon nanotube - is in contact with 
two electronic leads. Due to the tiny size of the struc¬ 
ture, single-electron tunneling processes can cause a tran¬ 
sient change of its electronic geometry. This change in 
combination with intermolecular interactions can cou¬ 
ple the electronic to the vibrational degrees of freedom. 
An important consequence is the appearance of non- 
linearities in the current-voltage characteristics and the 
conductanceThese effects can be associated 
with the possibility of inelastic processes due to the bias 
and to the formation of phonon sidebands in the excita¬ 
tion spectrum.— - — 

From a theoretical perspective, such a quantum 
dot setup can be described by the Anderson-Holstein 
mode l 19 i 20 . To its full extent, this model accounts for 
a tunneling coupling between the quantum dot and elec¬ 
tronic reservoirs, a linear coupling of the electrons oc¬ 
cupying the quantum dot to phonons, as well as an on¬ 
site Coulomb interaction. In this paper, we are mainly 
interested in the effects of the molecule’s vibration on 
the charge transport through the quantum dot so that 
it is expedient to consider a single vibrational mode and 
spinless electrons. Therefore, the model can be simplified 
to account for a single electronic level which is linearly 
coupled to a local phonon, whereas the electron-electron 
coupling is disregarded. In the framework of this spin¬ 
less Anderson-Holstein model, a lot of progress has been 
made, offering a deep insight into the physics caused by 
the coupling of the electron and the phonon (see, e.g., 
Refs. |15l - ll8ll2ll - |23j ). Besides approximative approaches 
numerically exact methods have (recently) become possi¬ 
ble for a nonequilibrium situation of a vibrational quan¬ 
tum dot (see, e.g., Refs. [23 - [28j ). 

A central quantity to describe the nonequilibrium 
transport through the quantum dot is its local density 


of states (LDOS). Single-particle observables such as the 
current or the dot occupation can be directly derived 
from il— . In equilibrium, the spectral density is well 
understood (see Refs. mmi and references therein). 
Moreover, the close connection between nanomechanical 
vibrations of the quantum dot and the sidebands has 
been confirmed. The interplay between nanomechanical 
vibrations and a finite bias voltage, however, still remains 
a challenging task outside certain limiting casesi 27 i 32 i 33 


In this paper, we address this problem in a numeri¬ 
cally exact way by using the diagrammatic Monte Carlo 
metho d 24 ' 34 ! 35 (diagMC). For this purpose, we use a two- 
terminal setup with an auxiliary electrode. This allows 
for an exact study of the LDOS on a quantum dot cou¬ 
pled to two electrodes in the limit of a vanishing coupling 
to the auxiliary lead— Throughout this paper we con¬ 
sider the deep quantum limit at T = 0. 


For a thorough discussion of the numerical results, 
we use an interpolative self-energy approximation (ISA), 
in which it is possible to include electron-phonon 
interaction s 18 ! 22 , and the well-established single parti¬ 
cle approximatio n 15 ! 16 ! 39- — (SPA). Although these ap¬ 
proximations rely on completely different approximation 
schemes, one uses the underlying common assumption 
that the phonons are described via an equilibrium dis¬ 
tribution. Consequently, effects due to a nonequilibrium 
distribution of phonons can be clearly identified by com¬ 
paring these approximations to the numerically exact re¬ 
sults. 


The structure of the paper is as follows: In SectionQ]we 
introduce the Anderson-Holstein model. In Section ITT1 we 
show how the diagMC can be used to calculate the LDOS 
for a quantum dot with an electron-phonon interaction 
by adapting the approach of Ref. [33. Moreover, the 
ISA, SPA and certain limiting cases are briefly summa¬ 
rized. The results for weak electron-phonon couplings are 
presented in Section eh The moderate polaronic regime 
is addressed in Section m 
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I. THE MODEL 

In our discussion, we consider a molecular quantum 
dot connected by a tunneling coupling to a left (L) and 
right (R) electrode. For such a system, a single-electron 
charging of the quantum dot can cause a nanomechan¬ 
ical vibration of the molecule^. A reasonable model 
Hamiltonian for this situation is provided by the spinless 
Anderson-Holstein mode l 42 ’ 43 (throughout this paper we 
use units with h = e = fcs = 1): 

H =J2 ( e «fc - do) °L°ak + (“1^ + dt «afc) 

fcGo: k£ot 

+ + Xd^d (b + &t) + ujob^b . (1) 

The quantum dot is modeled by a single electronic en¬ 
ergy level at £d- d) (d ) is the electron creation (annihi¬ 
lation) operator on the dot. a = L denotes the left and 
a = R the right electrode. The electronic creation (an¬ 
nihilation) operator on electrode a at energy level e a k is 
denoted by (a a k)- The bias voltage, V = /iL — Mr, 
is defined as the difference between the two chemical po¬ 
tentials n a of the respective electrode. This quantity is 
assumed to be constant for all times. 7 Q are the tunnel¬ 
ing amplitudes. The tunneling rates in the absence of 
manybody effects are given by T Q = 27r ( o Q |7 Q ,| 2 , where 
p a is the density of states of lead a, which is assumed 
to be a flat band. The linear coupling of the vibrational 
mode to the electronic degrees of freedom is described by 
a single phonon mode with frequency wo- b ^ and b are the 
phonon creation and annihilation operators, respectively. 
A denotes the electron-phonon coupling constant. 

In the subsequent discussion, we study the deep quan¬ 
tum limit at T = 0. Moreover, a symmetric setup is 
assumed where Tl = Tr = T/2, = —/tr = V/2 and 

we consider the particle-hole symmetric case, cd = 0. 
cd — cd — A 2 /wo is the polaron-shifted energy level of the 
quantum dot. This parameter regime is very interesting 
since the steady-state dot occupation is always (n) = 0.5 
due to the symmetry of the setup. Consequently, any 
bias dependence of the LDOS can only be caused by the 
electron-phonon interaction and not the electronic occu¬ 
pation of the quantum dot itself. 

II. APPROACHES FOR THE SPECTRAL 
DENSITY OF A VIBRATIONAL QUANTUM DOT 

A. Diagrammatic Monte Carlo simulation method 

Despite recent progress in developing analytical ap¬ 
proaches for the Anderson-Holstein model, e.g., by means 
of diagrammatic resummation schemes^—, a complete 
solution outside certain limiting cases is currently un¬ 
known. In order to calculate the spectral density with¬ 
out having to rely on methods which involve intrinsic ap¬ 
proximations, numerical methods are needed (see, e.g., 




FIG. 1: Sketch of the three-terminal setup with a weak cou¬ 
pling of the auxiliary lead to the quantum dot Tm -C IT, Tr. 
The bias voltage Vm can be tuned in order to access the de¬ 
sired energy in the spectral density of the quantum dot. 
is the respective chemical potential of the auxiliary electrode. 


Refs. [23 - [28j b A suitable approach to access regimes of 
arbitrary voltage, electron-phonon interaction, and dot- 
lead coupling strength is the numerically exact diagMC 
metho d 24 ’ 34 ’ 35 ’ 47 , which is able to simulate finite temper¬ 
atures as well as T = 0. 

In the subsequent discussion we use a similar approach 
to that of Ref. [H], where the diagMC has been used to 
calculate the LDOS for the Anderson impurity model. 
Here, we briefly summarize this approach and show how 
it can be adapted for the case of a local electron-phonon 
interaction on the quantum dot. 

Following the lines of Refs. [36H37j . the spectral den¬ 
sity of a two-lead quantum dot can be calculated exactly 
using an auxiliary lead at chemical potential /im with a 
vanishing dot-lead coupling (see Fig. Q] for a sketch of 
the setup). Their basic idea is to generalize the Meir- 
Wingreen formula for the current^ to a three-terminal 
setup to obtain 


— p + p M J dwp D (w) 

x [r/ M (w) - Tl/l(w) - r R /R(w)] , (2) 

where the subscript M denotes the auxiliary lead and 
Pd(w) is the LDOS of the quantum dot. / a (w) = 
l/( e /3(^-M Q ) _|_ j) denotes the Fermi function of lead a, 
where /3 is the inverse temperature. 

In the limit of a vanishing tunneling coupling of the 
auxiliary electrode, i.e. Tm —>• 0, one obtains the ther¬ 
mally broadened LDOS of the two-electrodes setup by 
deriving Eq. © with respect to the chemical potential of 
the auxiliary lead^i 
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lim 

Tm ->0 


- M 


9Im 

Ohm 


= J dwpD(w) 


<9/mM 

d/iM 


( 3 ) 


In the deep quantum limit, i.e. at T = 0, the derivative 
of the Fermi function in Eq. d3j) becomes a delta distribu¬ 
tion so that the exact LDOS of the two-electrodes setup 
is obtained and the thermal broadening vanishes^: 


P d(Mm)= lim (4) 

Fm —>0 0(1 m 

A convenient way to evaluate Eq. Q for the Anderson- 
Holstein model is to use the diagrammatic expansion in 
the tunnel couplin g 24 i 35 i 47 . This expansion allows for a 
complete decoupling of the influence of the leads to the 
dot, denoted by £m(s 7 j), from the phononic one includ¬ 
ing the dot’s energy level, denoted by Q(s n ). Therefore, 
with the use of Eq. © one obtains the transient which 
establishes the LDOS starting from an initially decoupled 
preparation: 


s n . For the considered single phonon mode, the bath 
autocorrelation function is given by 


L(s) = — [cos(wos) — isin(wos)] . (9) 

Wo 

We would like to emphasize that the only assumption 
for this approach is that the quantum dot is initially de¬ 
coupled from the leads. In detail, right before the cou¬ 
pling, the leads and the quantum dot are considered to 
be in their respective thermal equilibrium. The nonequi¬ 
librium aspect of the system enters via the coupling of 
the leads to the quantum dot at t = 0. This causes some 
transient dynamics until a nonequilibrium steady state is 
reached. 

Since Q(sn) is independent of the leads, the derivative 
with respect to the chemical potential of the auxiliary 
lead in Eq. 10 only acts on Cm(sn), which is a deter¬ 
minant of a matrix consisting of lesser and greater self¬ 
energies of the decoupled leads. Calculating this deriva¬ 
tive, one obtains at T = 0i2£ 


OO « 

Pd(mm) = 2 lim lim T^ 1 Y] (-1)" / ds* 
t-»oor M ->o J 


X -T — Re{£ M (Sn)<?(Sn)} ■ (5) 

OH M 


We used the abbreviation 


t 


0 


ds„ 



( 6 ) 


where s n = (si, s 2l ■ ■ ■, S 2 n} is the time-ordered sequence 
of 2 n tunneling times Sj. 

While in Ref. [38j Q accounts for a Coulomb on-site 
interaction, in Eq. © it provides the influence of the 
electron-phonon interaction given by21 


g(s n ) = T[s n ] ^{s,-s 2 +s 3 --) ) ( 7 ) 

where an initially empty quantum dot is consid¬ 
ered. ~F[s ra ] denotes the Feynman-Vernon influence 
functional^: 


J"[s„] = exp 



>C:s2<s\ 


ds 2 g(si)L(si - s 2 )q(s 2 ) 


( 8 ) 


The integrations are performed on the Keldysh contour 
C : 0 —^ t —0. q(s) denotes the occupation of the 
quantum dot at time s which is fully determined by the 
initial condition of the quantum dot and the position 
as well as the number of the tunneling events given by 


lim 

r M ->o Fm o/im 


Cm (s„) = i” det(<S M (sn)) 


( 10 ) 


where 


{ £ < (s 2 fc-i,s 2:) ) , for j < k 

£ > (s 2 fc-i,s 2 j) , for j > k 

_L, if V s 2j ) = t 

( 11 ) 

We would like to emphasize that the limit Tm —> 0 
has been performed analytically so that the auxiliary 
lead is not influencing the two-terminal quantum dot. 
Consequently, the stationary limit of Eq. © is the ex¬ 
act LDOS of the two-terminal setup. We would like to 
note that the steady state is defined with respect to the 
reduced dynamics of the quantum dot. Therefore, the 
complete system is in nonequilibrium even though a time- 
independent steady state for observables on the quantum 
dot is reached. 

Therefore, for a given sequence of tunneling events, s n , 
it is straightforward to calculate the influence of the leads 
on the dot as well as the phononic influence without any 
approximation. The summation and integration over all 
possible tunneling events in Eq. © can be done con¬ 
veniently in a numerical exact manner by using Monte 
Carlo samplin g 24 i 34 i 35 . Using this method, the only oc¬ 
curring error is a controllable statistical one. We note 
that whenever no error bar is visible in the subsequent 
figures, the error is smaller than the symbol size. 

In the following we will consider two different coupling 
procedures of the leads to the dot at t = 0: A sudden, and 
a smooth switch-on (for details of the coupling procedure 
see, e.g., Ref. Hi). In addition, we truncate the leads’ 
density of states at a value ±e c . The reason for this sharp 
cutoff is that an instantaneous coupling of the electrodes 
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to the quantum dot can lead to excitations in the leads, 
which are arbitrarily high in energy at t = 0^2 These 
short-living excitations are not only unphysical but also 
make a numerical evaluation using diagMC unfeasible. 
For our results, the cutoff is chosen to be the largest 
energy scale in the system so that a further increase of 
e c does not change our results for times t > 


B. Limiting cases 


In this section, we briefly discuss two limiting cases, 
which can be solved analytically. The first one is the 
absence of electron-phonon interactions, A/T —> 0. Here, 
it is straightforward to see that the LDOS is independent 
of the applied bias voltage^: 


Pd(w) 


i r 

2^ (w - e D ) 2 + (T/2 f ' 


( 12 ) 


Considering a weakly coupled phonon mode, the electron- 
phonon coupling can be treated perturbativel y 16 1 21 ■ 51 Jil, 
Consequently, Eq. indicates that in the perturba¬ 
tive regime no, or only a weak voltage dependence in the 
spectral density is expected. 

Similar arguments hold for the atomic limit, where the 
electron-phonon coupling becomes very large, T/A —> 0. 
In this case, the LDOS is given by sharp delta-peaks at 
multiples of the phonon frequency^ 


OO fc 

Pv(uj) = e~ 9 ^2 it 1 - («d)] S(iJ - e D - ku 0 ) 
k =0 

+ (n D )<5(w - e D + kuj 0 )} , (13) 

where (tt-d) is the charge on the quantum dot, and 
g = (A/wo) 2 is the dimensionless electron-phonon cou¬ 
pling strength. Since in this limit the leads are decoupled 
from the quantum dot, the LDOS cannot depend on the 
bias voltage. 

To summarize these considerations, it is clear that in 
the weak as well as in the strong coupling limit, a possi¬ 
ble voltage dependence of the LDOS can only be weak. 
Therefore, the nonequilibrium LDOS outside these lim¬ 
iting cases is expected to be the most interesting one. 


C. Approximative approaches 


Finally, we discuss two important and often used ap¬ 
proximative approaches. 

A very popular approximate scheme is the 
S P A 15 i 16 i 31 i 39 — , where the electrons are decoupled 
from the phonons. For this approach the LDOS at 
temperature T can be written as 


F °° 

= —e-« (2 ” B+1) £ 4[2 5 V(n B + l)n B ]e fe ^°/ 2 

k=—oo 


1 - (n D ) 


(to - e d — kcoo) 2 + (T/2) 2 


+ - 


W) 


(w — + kujo) + (r/ 2 ) 


(14) 


where the charge on the quantum dot, (nu), has to 
be calculated self-consistently. 4 is the modified Bessel 
function, /3 = 1/T and n B = l/(e^ w ° — 1). The simple 
structure of the SPA allows for a straightforward eval¬ 
uation. Moreover, it provides good results if the cor¬ 
relations between electrons and phonons are weak, or 
if the quantum dot is either empty or occupied. Fur¬ 
thermore, the atomic limit as well as the case of absent 
phonons is recovered. Outside these limiting cases, meth¬ 
ods are needed which go beyond the simple SPA decou¬ 
pling scheme. A well-established method, which is known 
to provide reasonable results for a broad range of param¬ 
eters is the ISA 22 . The basic idea is to perform a func¬ 
tional interpolation of the self-energies from the weak 
to the strong coupling regime. This scheme was origi¬ 
nally derived for the Anderson impurity mode l 56 ’ 57 and 
has been extended and widely used in different systems: 
multilevel quantum dots^, out-of-equilibrium transport 
through a single leve l 59 ’ 60 , and in dynamical mean-field 
theor y 61 ’ 62 to analyze the Mott transition in Hubbard- 
like models. For the nonequilibrium Anderson-Holstein 
model this approach provides accurate results beyond 
perturbation theory or SPAi£. 

Using the ISA and the SPA, the effects of the electron- 
phonon interaction on the charge transport can be dis¬ 
cussed qualitatively as long as their basic underlying as¬ 
sumption is fulfilled: an equilibrium distribution of the 
phonons. Since both methods cover a broad range of pa¬ 
rameters this in turn implies that if the physics is quali¬ 
tatively not covered by either of these methods, there is 
a strong indication that the phonons no longer obey an 
equilibrium distribution. 


III. WEAKLY COUPLED PHONON MODE 

We start the discussion of our results by considering 
the weak-coupling regime with A = T, and a rather large 
phonon frequency wo = 4T. The corresponding polaronic 
self-energy can be determined to be A po i = A 2 /wo = F/4. 
Since A po i < T the formation of a polaron is a relatively 
rare event and the electrons are thus weakly coupled to 
the phonons. 

In Fig. [2] the transients which establish the spectral 
density for V = 2T at to = ±0.5r are shown. An instan¬ 
taneous coupling to the electrodes leads in this case to 
an overshooting, with the steady state being approached 
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FIG. 2: Time-dependent diagMC results for po(u},t) with 
A = F, tuo = 4T for V ~ 2T and lo = 0.5r (red circles) as 
well as to = — 0.5T (blue diamonds). Empty symbols denote 
results for a switch-on time of r sw = lOE^ 1 . The bandwidth 
is 2e c = 8T. 
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monotonically. The relevant timescales for the dynam¬ 
ics can be estimated to be O (T _1 ). A smooth coupling 
of the quantum dot to the electrodes^ establishes the 
steady state adiabatically. 

Since the particle-hole symmetric case is considered, 
the spectral density in the stationary limit must be sym¬ 
metric with respect to w = 0. In Fig. [2] it can be seen 
that in the weak-coupling case this property is well ful¬ 
filled even in the transient regime. 

The resulting LDOS in the frequency domain is shown 
in Fig. [3] for two different voltages: V = 2T and V = 
10T. Comparing results from ISA to the results extracted 
from the time-dependent diagMC calculations, we find 
excellent agreement. 

The overall shape of the spectral density is very sim¬ 
ilar for both voltages. A comparison with the results in 
the absence of phonons (A = 0) reveals that the height 
of the central peak remains almost unchanged in the 
low-voltage regime. The ISA reveals a slight decrease 
of the central peak when increasing the bias voltage. 
Since this decrease is small, the Friedel-Langreth sum 
rul o 63 ^ 64 , which pins the height of the central peak to 
Pd(0) = 2/('kT)^ provides a good approximation also 
for the nonequilibrium situation. 

Small sidebands at multiples of the phonon frequency 
are observed, which are independent of the voltage within 
the accuracy of the results extracted from diagMC. The 
ISA reveals features at |w| = fcwo ± V/2, with k being an 
integer, where the LDOS is changing rapidly (see inset 
of Fig. [3]). These features can be attributed to inelas¬ 
tic electron tunneling processes as predicted by different 
theoretical approache s 18 ! 44 to appear both in the spectral 
density and conductance. For the considered nonequi- 
librium spectral density with a particle-hole symmetric 
setup they appear at the condition |w| = fcwo ± V/2. For 
large k the effect can be tiny due to the large amount of 
phonons involved. 

We conclude that the voltage dependence of the spec¬ 
tral density for a weakly coupled phonon mode is very 
small. This observation is in excellent agreement with 
our discussion of the limiting cases given in Section III Bl 



<oin 


FIG. 3: LDOS at zero temperature for a weakly coupled 
phonon mode with A = T, uo = 4F. Red indicates V = 2F 
and green V = 10T. The straight lines are the ISA, circles 
and diamonds denote the diagMC results. The LDOS in the 
absence of phonons (A = 0), given by Eq. (1121) . is shown as a 
black dashed line. The inset shows a zoom into the region of 
the first phonon sideband at w = wo = 4T. 


An important consequence for future theoretical ap¬ 
proaches is that in this regime it is sufficient to solve the 
nonequilibrium problem by calculating the LDOS using 
equilibrium theory. The nonequilibrium aspect of the 
system enters only via the integration limits for single¬ 
particle observables such as the dot occupation or the 
current. 


IV. SPECTRAL DENSITY IN THE MODERATE 
POLARONIC REGIME 

In the polaronic regime, the formation time of a po- 
laron is shorter than the average occupation time of the 
electron on the quantum dot. Consequently, the for¬ 
mation of a polaron becomes likely so that pronounced 
phonon sidebands are expected. The corresponding pa¬ 
rameter regime can be determined to be A po i > F. 

The focus of the subsequent discussion is on the volt¬ 
age dependence of the spectral density. Our approach 
to distinguish different voltage regimes is by considering 
the number of phonon sidebands, which are included in 
the voltage window: If only the central transport channel 
is between the two chemical potentials, that is V < ujo, 
the low-voltage regime is realized. For voltage windows 
including one or more sidebands, we expect that nonequi- 
librium aspects are most pronounced. 


A. Low-voltage regime 

Important references in the low-voltage regime are 
equilibrium results such as the Friedel-Langreth sum 
rule, which is fulfilled at V = 0 independent of the 
electron-phonon interaction strength—^. Since phonon 
sidebands form when increasing the electron-phonon 
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FIG. 4: Time-dependent diagMC results for the central reso¬ 
nance at w = 0 with V = 2T for various phonon parameters: 
A = 3T, u>o = 4r (yellow triangles), A = 2F, wo = 2F (red 
circles), A = 4F, u>o = 3T (green diamonds), and A = 8F, 
uj o = 4F (blue pentagons). The bandwidth is set to 2e c = 8T. 
The corresponding transients in the spirit of Ref. [46l| . given 
by Eq. HU), are lines with the same color code as the respec¬ 
tive numerical data. The analytical result of the transient in 
the absence of phonons is shown as a black line. 
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FIG. 5: diagMC results for pT>(uo,t) for A = wo = V = 2F. 
The upper panel shows w = 0.5F, the lower one w = F. Empty 
circles denote a smooth switching within r sw = 12r _1 , and 
filled ones correspond to an instantaneous coupling at t = 0. 
The bandwidth is 2e c = 8T. 


couplin g 16 i 20 i 30 i 31 , and the norm of the spectral density 
needs to be preserved, the central peak must be phonon- 
narrowed. Such a narrowing can be determined to be 
r ~ re -ff £° This narrowing of the central peak for strong 
electron-phonon couplings is the origin of the Franck- 
Condon blockade effect discussed in detail in Ref. 65j 
Regarding the nonequilibrium problem, in Ref. 46( it 
was shown via an approximative study that such a nar¬ 
rowing of the resonances can increase the timescales rel¬ 
evant for the charge transport up to 


T pol = exp 


(X/uoY 


T -i = f-i i 


(15) 


Correspondingly, e.g., the transients for establishing the 
central peak of the spectral density at uj = 0 can be 
estimated to be: 

p D (w = 0, t) = ^ (l - e“ tr/2 ) . (16) 

Our diagMC results confirm this behavior as shown 
in Fig. U for various electron-phonon couplings in the 
moderate polaronic regime. We would like to empha¬ 
size the broad range of phonon parameters, for which 
Eq. (1161) provides accurate results. Besides confirming 
the existence of phonon-induced long timescales^ deter¬ 
mined by r _1 , an important consequence of this behavior 
is that in the low-bias regime the Friedel-Langreth sum 
rule is fulfilled. That is, the central peak is pinned to 
Pd(w = 0) = 2/(7rr) independent of the electron-phonon 
coupling strength. 

In the subsequent discussion, we will study a param¬ 
eter regime where nonequilibrium effects are most pro¬ 
nounced. According to our preceding discussion, it must 
therefore be neither close to the limiting cases of a very 
weak phonon coupling nor in the strong polaronic regime. 


Another requirement is that we are able to extract the 
complete LDOS from the time-dependent diagMC re¬ 
sults. Therefore, the steady state has to be reached 
within times which are accessible by diagMC. A reason¬ 
able choice of parameters fulfilling these requirements is 
A = wo = 2T: the moderate polaronic regime is accessed 
and the longest timescales of the transients are roughly 
T” 1 ft: 2.7T~ 1 . The current implementation of the di¬ 
agMC is able to simulate up to t ~ C>(10r -1 ) within 
reasonable computational effort. Consequently, we can 
not only discuss the transient dynamics, but it is also 
reasonable to extract the steady state of the system from 
the time-dependent results. 

In Fig. [5] the transients which establish the spectral 
density at u = 0.5T and ui = T are shown. Similar to the 
weak-coupling regime, an instantaneous coupling of the 
electrodes to the quantum dot leads to an overshooting 
for small times. The steady state, however, is approached 
non-monotonically in an oscillating manner. The charac¬ 
teristic timescale for the convergence towards the steady 
state can be estimated to be given by F _1 , pronounced 
features in the oscillations occur with a periodicity given 
by t ~ 27t/wo- These transients can be reduced by a 
smooth switch-on procedure of the leads to the quantum 
dot. If a rather long switch-on time is used, the steady 
state can be extracted with good accuracy. 

It is worth noticing an interesting effect that is clearly 
observed in the transient regime for the strong-coupling 
case: while the particle-hole symmetry requires that 
the stationary spectral density is symmetric, pd(w) = 
Pd(— w), this relation doesn’t have to be fulfilled neces¬ 
sarily in the transients. In Fig. [6] an overshooting for 
a positive frequency, uj = ujq, can be observed, whereas 
the transients for a negative frequency, u> = — ujo, show a 
monotonic increase in time. This effect can be explained 
by the asymmetry in the initial preparation: right be¬ 
fore the coupling of the quantum dot to the leads, the 















FIG. 6: Same plot as in Fig. 0 however, for the first phonon 
sideband with to = too = 2Y (red circles) and to = —too = 
—2r (blue diamonds). The green triangles denote the diagMC 
results of the average pD,av(w,t), as in Eq. (11711 . A smooth 
switching within r aw = 6W 1 is employed. 


quantum dot is empty. Therefore, only resonances posi¬ 
tive in energy exist, since no deexcitation of phonons is 
possible at t = 0. On the finite timescale necessary to 
establish the spectral density^ this asymmetry in the ini¬ 
tial preparation leads to an asymmetry in the transients. 
Combined with the phonon-induced long timescales, this 
fact provides a deeper understanding of the splitting of 
the current depending on the initial preparation on a 
timescale given by O (r _1 ), which was observed, e.g., in 
Ref. [66} . This is further corroborated by the fact that for 
the considered particle-hole symmetric case one obtains 
p^n pty (a;,t) = p° ccupied (_ w , f), where the superscript de¬ 
notes the initial occupation of the quantum dot. 

Since neither po(ui = u>o,t), nor pd(w = —u>o,t) ex¬ 
hibit a clear steady state in Fig. [Gl we use the particle- 
hole symmetry of the considered setup, and define the 
average of the transients by 

/9D,av(w,i) = ^ [pn(u,t) +p D (-u>,t)] • (17) 

The particle-hole symmetry ensures that this func¬ 
tion has the same steady-state value as pn(w,t) and 
Pd(— w,t), separately. The transient dynamics, however, 
show a quicker convergence towards the plateau value due 
to the averaging between excitations and deexciations in 
Eq. HT|. Therefore, the steady state of pD,av(w,t) can 
be extracted with reasonable accuracy. We would like 
to note that similar convergences of the observables can 
also be found for the currents: While the average current 
reaches a plateau for t > 8T _1 , the left and the right cur¬ 
rent converge to a joint steady state for times t > lir -1 . 

In Fig. m we plot the extracted spectral density of 
the quantum dot, and we make a comparison between 
diagMC, ISA, as well as SPA. It should be stressed that 
the error bar of the extracted steady state from the time- 
dependent diagMC results is twice the total change of 
dD,av(w,t) from t = 8r -1 to t = 12r _1 . 

A very sharp central peak at to = 0 is observed with 
a height given by the Friedel-Langreth sum rule. Com¬ 
pared to the case of absent phonons, the width of the 
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FIG. 7: Zero temperature spectral density for A = too = 2F 
and V = 2r. The results extracted from diagMC are shown 
with red filled circles, the ISA is a blue line, and SPA is green. 


central peak is reduced to f s=a Te~ 9 , as it was also ob¬ 
served in Ref. [3(J for the equilibrium situation. This 
observation eventually proves the close connection be¬ 
tween the long transients and phonon narrowing of the 
resonances, which was discussed in Ref. [46j] by means of 
an approximative method. 

Phonon sidebands can be found at multiples of the 
phonon frequency, with an exponentially decreasing 
height. This behavior reflects the fact that transport 
outside the voltage window is strongly suppressed. More¬ 
over, clear dips in the spectral density appear between 
two phonon sidebands. 

The ISA describes the results obtained from diagMC 
with remarkably good accuracy. Small differences are 
only visible at the first phonon sidebands, where the ISA 
predicts a slightly larger value. Since one basic assump¬ 
tion of the ISA is an equilibrium phonon distribution, 
we conclude that effects due to a (possible) nonequilib¬ 
rium phonon distribution only play a minor role in the 
moderate polaronic regime at low biases. 

Regarding the low-frequency domain calculations us¬ 
ing the SPA a clear deviation from the diagMC is ob¬ 
served. This means that the charge on the quantum dot is 
strongly correlated with the excitation of phonons for fre¬ 
quencies that are not too large. In the large frequency do¬ 
main, however, a good agreement between the diagMC, 
and the SPA is observed. Therefore, the electrons are 
decoupled from the phonons for transport with energies 
much larger than the voltage window, which confirms the 
results of Ref. HH . 

B. Far-from-equilibrium spectral density 

In the subsequent discussion, we will analyze nonequi¬ 
librium effects in the moderate polaronic regime. For this 
purpose we consider voltage windows that contain one or 
more phonon sidebands. While nonequilibrium effects 
are not important in the weak-coupling regime, as dis- 
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FIG. 8: Time-dependent diagMC results for the transients, 
which establish the central peak pd(w = 0) for A = wo = 2F 
and various voltages: V = 2F (red circles, bandwidth 2e c = 
8F), V = 6r (green triangles, 2e c = 8T), V = 10T (blue 
diamonds, 2e c = 12T), V = 14T (yellow reverted triangles, 
2e c = 16r). The approximative description in the spirit of 
Ref. [4(|, given by Eq. (fTSll . is shown as a red line. 


cussed in Section cm the effect of a large bias voltage is 
strong in the moderate polaronic regime: The transients, 
which establish the central transport channel at ui = 0, 
are shown in Fig. [HI where the same phonon parameters 
are used as in the previous section: A = ojq = 2T. For 
a voltage window, which includes one or more phonon 
sidebands, the transients no longer follow the exponential 
convergence, given by Eq. m- Moreover, the relevant 
timescales for the transient dynamics are significantly 
smaller so that the steady state is reached faster than 
in the low-bias regime. Furthermore, the steady-state 
value drops to a much smaller value, which reflects the 
fact that the transport through the quantum dot outside 
the low-voltage regime is dominated by the excitation of 
one or more phonons. This behavior clearly violates the 
Friedel-Langreth sum rule. 

The strongest voltage dependence of the height of the 
spectral density at w = 0 is observed when increasing 
the voltage from V = 2F to V = 6r. The reason for this 
behavior is that for A = wo = 2T the central transport 
channel as well as the first phonon sideband are most pro¬ 
nounced for V = 2T as it can be seen in Fig. [7] Includ¬ 
ing the first sideband into the voltage window by setting 
V = 6T, charge transport involving single-phonon pro¬ 
cesses becomes likely and thus this important transport 
channel opens. This reduces the probability for charge 
transport without exciting a phonon and, consequently, 
the height of the central resonance at w = 0 decreases. A 
further increase of the voltage does not exhibit this pro¬ 
nounced behavior since the transport channels involving 
two or more phonons are much smaller for the considered 
parameters. 

An important consequence of the decreasing central 
transport channel is that the weight of the LDOS is 
shifted towards larger frequencies. A similar shift has 
been reported recently in the differential conductance^. 
This shift causes an increase of the phonon sidebands 



FIG. 9: Time-dependent diagMC results for pD,av(w,t) with 
|ta| = wo (upper panel) and |w| = 3wo (lower panel) with A = 
wo = 2r. V = 2T is highlighted with red circles (bandwidth 
2e c = 8r in the upper panel, and 2e c = 14F in the lower one), 
and V = 10r with blue diamonds (bandwidth 2e c = 14T 
for both panels). A smooth switching within r sw = 6F _1 is 
employed. 


outside the voltage window as it can be seen in the lower 
panel of Fig. [9] Including a phonon sideband into the 
voltage window causes a drop of the phonon resonance 
as it can be seen in the upper panel of Fig. [9] 

The resulting spectral density extracted from the time- 
dependent diagMC for V = 10T is shown in Fig. [Tol 
A comparison with the low-bias LDOS in Fig. [3 reveals 
that inside the voltage window all peaks seem to align 
to a similar height, whereas peaks outside the voltage 
window are increased. Moreover, the narrowing of the 
central transport channel to T ~ e _s r is no longer ob¬ 
served in the large voltage regime. Rather, a width of 
approximately T, which is the value for the interaction- 
free case, is recovered. We note, however, that this is 
not an indication that the charge transport through the 
quantum dot is uncorrelated from the phonons: a com¬ 
parison with the SPA reveals that the central transport 
channel is strongly suppressed. 

In addition, Fig. [lO] reveals a clear deviation of the 
ISA from the diagMC results. For such a large voltage 
the ISA spectral density has (almost) converged to the 
SPA case. It is interesting to remark that a recent di¬ 
agrammatic resummation scheme valid in the polaronic 
regime^ predicts a similar convergence towards the SPA 
for large voltages. As these approximate theories do not 
include the effect of the nonequilibrium distribution of 
phonons in a self-consistent way, the numerically exact 
diagMC results strongly indicate that this effect is impor¬ 
tant in the moderate polaronic regime in the large bias 
limit. 

A straightforward way to confirm that the deviations 
observed in Fig. [TU] are indeed produced by a nonequilib¬ 
rium phonon distribution is to try to simulate the results 
by an equilibrium distribution at some effective temper¬ 
ature. For this purpose, we use the SPA since the ISA as 
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FIG. 10: Same color code as Fig. [7] but with V = 10F. 
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FIG. 11: Same diagMC results as in Fig. [TO] The SPA has 
been calculated with an effective temperatures of the phonons: 
T = F (blue), T = 3F (black), and T = 10F (green). 


well as the approach of Ref. 0 converge to the SPA for 
sufficiently large bias voltages. 

In Fig.[lT]the SPA results given by Eq. (TTH) for various 
effective phonon temperatures T —> T e g are compared to 
the diagMC results, which are calculated for V = 10r 
and T = 0. Strikingly, the central peak of the SPA spec¬ 
tral density decreases by increasing the phonon temper¬ 
ature, whereas phonon resonances well away from the 
central peak increase - a behavior similar to the diagMC 
results for increasing voltage. For an effective phonon 
temperature of roughly T e g ss 3r, the results from SPA 
match the diagMC results for T = 0 and V = 10r with 
good accuracy. We would like to note that despite the 
good overall agreement, small differences can be observed 
for the second phonon resonance, which means that not 
all effects can be completely described in detail by this 
effective theory. 


V. CONCLUSIONS 

In this paper we calculated the nonequilibrium spectral 
density of a vibrational quantum dot using the numerical 
exact diagMC technique and compared its predictions 
with those of approximate methods such as ISA and SPA. 

We showed that for a weak electron-phonon interaction 
the spectral density of the quantum dot resembles the 
equilibrium one independently of the bias voltage. For 
intermediate electron-phonon coupling strengths in the 
moderate polaronic regime, we determined a significant 
voltage dependence of the spectral density. An increas¬ 
ing bias voltage shifts the weight of the spectral den¬ 
sity towards larger energies: The central resonance de¬ 
creases, whereas phonon resonances outside the voltage 
window are increased with respect to the equilibrium re¬ 
sults. Inside the voltage window our results indicate that 
the phonon peaks align to a similar height. We were able 
to link the voltage dependence of the spectral density to 
an effective “heating” of the phonons caused by inelastic 
excitations. 

The explicit voltage dependence of the spectral den¬ 
sity points out the importance of accessing the spec¬ 
tral density directly, e.g., by means of a three-terminal 
setup 36 i 37 i 67 . An indirect measurement, e.g., of the differ¬ 
ential conductance might lead to a discrepancy between 
the result and the actual spectral density due to its volt¬ 
age dependence. 

Another consequence of our findings is that for future 
descriptions by means of approximative approaches it is 
desirable to also account for nonequilibrium effects of the 
phonon distribution, which could be preformed, e.g., as 
proposed in Ref. |32| . 

Finally, we would like to emphasize that for small volt¬ 
ages, we confirmed the existence of phonon-induced long 
transients previously proposed in Ref. [46j. Moreover, 
we pointed out that the inverse of the width of the res¬ 
onances in the spectral density determines the relevant 
timescale in the system. 

Acknowledgments 

The authors like to thank R. C. Monreal, A. Levy Yey- 
ati, R. Seoane Souto, and A. Komnik for many fruitful 
discussions. KFA acknowledges the computing time at 
the bwGRID and Juropa in Jiilich. This work was fi¬ 
nancially supported by Spanish Mineco through grant 
FIS2011-26516, and by the NSF (PIF-1211914 and PFC- 
1125844). 


1 H. Park, J. Park, A. K. L. Lim, E. H. Anderson, A. P. 
Alivisatos, and P. L. McEuen, Nature 407, 57 (2000). 

2 R. H. M. Smit, Y. Noat, C. Untiedt, N. D. Lang, M. C. 


van Hemert, and J. M. van Ruitenbeek, Nature 419, 906 

( 2002 ). 

3 N. B. Zhitenev, H. Meng, and Z. Bao, Phys. Rev. Lett. 88, 









10 


226801 (2002). 

4 X. H. Qiu, G. V. Nazin, and W. Ho, Phys. Rev. Lett. 92, 
206102 (2004). 

5 J. Kushmerick, J. Lazorcik, C. Patterson, R. Shashidhar, 

D. S. Seferos, and G. C. Bazan, Nano Lett. 4, 639 (2004). 

6 N. Liu, N. Pradhan, and W. Ho, J. Chem. Phys. 120, 
11371 (2004). 

7 L. H. Yu, Z. K. Keane, J. W. Ciszek, L. Cheng, M. P. 
Stewart, J. M. Tour, and D. Natelson, Phys. Rev. Lett. 
93, 266802 (2004). 

8 A. N. Pasupathy, J. Park, C. Chang, A. V. Soldatov, 
S. Lebedkin, R. C. Bialczak, J. E. Grose, L. A. K. Donev, 

J. P. Sethna, D. C. Ralph, et al., Nano Lett. 5, 203 (2005). 

9 S. Sapmaz, P. Jarillo-Herrero, Y. M. Blanter, C. Dekker, 
and H. S. van der Zant, Phys. Rev. Lett. 96, 026801 (2006). 

10 R. Leturcq, C. Stampfer, K. Inderbitzin, L. Durrer, C. Hi- 
erold, E. Mariani, M. G. Schultz4, F. von Oppen, and 

K. Ensslin, Nature Physics 5, 327 (2009). 

11 D. Seeker, S. Wagner, S. Ballmann, R. Hartle, M. Thoss, 
and H. B. Weber, Phys. Rev. Lett. 106, 136807 (2011). 

12 J. C. Cuevas and E. Scheer, Molecular Electronics: An 
Introduction to Theory and Experiment (World Scientific, 
Singapore, 2010). 

13 S. Ballmann, W. Hieringer, D. Seeker, Q. Zheng, J. A. 
Gladysz, A. Gorling, and H. B. Weber, Chem. Phys. Chem. 
11, 2256 (2010). 

14 J. Park, A. N. Pasupathy, J. I. Goldsmith, C. Chang, 
Y. Yaish, J. R. Petta, M. Rinkoski, J. P. Sethna, H. D. 
Abruna, P. L. McEuen, et al., Nature 417, 722 (2002). 

15 A.-P. Jauho, N. S. Wingreen, and Y. Meir, Phys. Rev. B 
50, 5528 (1994). 

16 K. Flensberg, Phys. Rev. B 68, 205323 (2003). 

17 A. Mitra, I. Aleiner, and A. J. Millis, Phys. Rev. B 69, 
245302 (2004). 

18 R. C. Monreal, F. Flores, and A. Martfn-Rodero, Phys. 
Rev. B 82, 235412 (2010). 

19 T. Holstein, Annals of Physics 8, 343 (1959). 

20 A. C. Hewson and D. Meyer, Journal of Physics: Con¬ 
densed Matter 14, 427 (2002). 

21 M. Galperin, A. Nitzan, and M. A. Ratner, Phys. Rev. B 
74, 075326 (2006). 

22 A. Martfn-Rodero, A. Levy Yeyati, F. Flores, and R. C. 
Monreal, Phys. Rev. B 78, 235112 (2008). 

23 T. Koch, J. Loos, A. Alvermann, and H. Fehske, Phys. 
Rev. B 84, 125131 (2011). 

24 L. Miihlbacher and E. Rabani, Phys. Rev. Lett. 100, 
176403 (2008). 

25 H. Wang and M. Thoss, J. Chem. Phys. 131, 024114 
(2009). 

26 R. Hiitzen, S. Weiss, M. Thorwart, and R. Egger, Phys. 
Rev. B 85, 121408 (2012). 

27 A. Jovchev and F. B. Anders, Phys. Rev. B 87, 195112 
(2013). 

28 L. Simine and D. Segal, JCP 138, 214111 (2013). 

29 Y. Meir and N. S. Wingreen, Phys. Rev. Lett. 68, 2512 
(1992). 

30 A. C. Hewson and D. M. Newns, Journal of Physics C: 
Solid State Physics 13, 4477 (1980). 

31 R. C. Monreal and A. Martfn-Rodero, Phys. Rev. B 79, 
115140 (2009). 

32 M. Galperin, A. Nitzan, and M. Ratner, Phys. Rev. B 73, 
045314 (2006). 

33 M. A. Laakso, D. M. Kennes, S. G. Jakobs, and V. Meden, 
New Journal of Physics 16, 023007 (2014). 


M. Schiro and M. Fabrizio, Phys. Rev. B 79, 153302 
(2009). 

E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, 
M. Troyer, and P. Werner, Rev. Mod. Phys. 83, 349 (2011). 

Q. -f. Sun and H. Guo, Phys. Rev. B 64, 153306 (2001). 

E. Lebanon and A. Schiller, Phys. Rev. B 65, 035308 

( 2001 ). 

L. Miihlbacher, D. F. Urban, and A. Komnik, Phys. Rev. 
B 83, 075107 (2011). 

D. M.-T. Kuo and Y. C. Chang, Phys. Rev. B 66, 085311 

( 2002 ). 

U. Lundin and R. H. McKenzie, Phys. Rev. B 66, 075303 

( 2002 ). 

A. Zazunov and T. Martin, Phys. Rev. B 76, 033417 
(2007). 

L. Glazman and R. Shekhter, Sov. Phys. JETP 67, 163 
(1988). 

M. Galperin, M. A. Ratner, and A. Nitzan, J. Phys.: Con- 
dens. Matter 19, 103201 (2007). 

R. Seoane Souto, A. L. Yeyati, A. Martfn-Rodero, and 

R. C. Monreal, Phys. Rev. B 89, 085412 (2014). 

S. Maier, T. L. Schmidt, and A. Komnik, Phys. Rev. B 83, 
085401 (2011). 

K. F. Albrecht, A. Martfn-Rodero, R. C. Monreal, 

L. Miihlbacher, and A. Levy Yeyati, Phys. Rev. B 87, 
085127 (2013). 

P. Werner, T. Oka, and A. J. Millis, Phys. Rev. B 79, 
035320 (2009). 

R. P. Feynman and F. L. Vernon, Jr., Ann. Phys. (Leipzig) 
24, 118 (1963). 

T. L. Schmidt, P. Werner, L. Miihlbacher, and A. Komnik, 
Phys. Rev. B 78, 235110 (2008). 

P. W. Anderson, Phys. Rev. 124, 41 (1961). 

F. Haupt, T. Novotny, and W. Belzig, Phys. Rev. Lett. 
103, 136601 (2009). 

T. L. Schmidt and A. Komnik, Phys. Rev. B 80, 041307 
(2009). 

R. Avriller and A. Levy Yeyati, Phys. Rev. B 80, 041309 
(2009). 

R.-P. Riwar and T. L. Schmidt, Phys. Rev. B 80, 125109 
(2009). 

G. Mahan, Many-particle physics (Plenum press, New 
York, 1993). 

A. Martin-Rodero, F. Flores, M.Baldo, and R. Pucci, Solid 
State Commun 44, 911 (1982). 

A. Martin-Rodero, E. Louis, F. Flores, and C. Tejedor, 
Phys. Rev. B 33, 1814 (1986). 

A. L. Yeyati, F. Flores, and A. Martin-Rodero, Phys. Rev. 
Lett 83, 600 (1999). 

A. L. Yeyati, A. Martin-Rodero, and F. Flores, Phys. Rev. 
Lett 71, 2991 (1993). 

A. A. Aligia, Phys. Rev. B 74, 155125 (2006). 

H. Kajueter and G. Kotliar, Phys. Rev. Lett 77, 131 
(1996). 

M. Pothoff, T. Wegner, and W. Nolting, Phys. Rev. B 55, 
16132 (1997). 

J. Friedel, Philosophical Magazine 7, 153 (1951). 

D. C. Langreth, Phys. Rev. 150, 516 (1966). 

J. Koch, F. von Oppen, and A. V. Andreev, Phys. Rev. B 
74, 205438 (2006). 

K. F. Albrecht, H. Wang, L. Miihlbacher, M. Thoss, and 
A. Komnik, Phys. Rev. B 86, 081412(R) (2012). 

R. Leturcq, L. Schmid, K. Ensslin, Y. Meir, D. C. Driscoll, 
and A. C. Gossard, Phys. Rev. Lett. 95, 126603 (2005). 


34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 

61 

62 

63 

64 

65 

66 

67 


